// ///////////////////////////// MIT License //////////////////////////////////// //
//                                                                                //
// Copyright (c) 2010 David Zsolt Manrique                                        //
//                    david.zsolt.manrique@gmail.com                              //
//                                                                                //
// Permission is hereby granted, free of charge, to any person obtaining a copy   //
// of this software and associated documentation files (the "Software"), to deal  //
// in the Software without restriction, including without limitation the rights   //
// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell      //
// copies of the Software, and to permit persons to whom the Software is          //
// furnished to do so, subject to the following conditions:                       //
//                                                                                //
// The above copyright notice and this permission notice shall be included in     //
// all copies or substantial portions of the Software.                            //
//                                                                                //
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR     //
// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,       //
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE    //
// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER         //
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,  //
// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN      //
// THE SOFTWARE.                                                                  //
//                                                                                //
// ////////////////////////////////////////////////////////////////////////////// //

//!
//! moltool
//!

// standard headers
#include <iostream>
#include <sstream>
#include <fstream>
#include <string>
#include <vector>
#include <list>
#include <set>
#include <map>
#include <iomanip>
#include <cmath>
#include <cstdarg>
#include <stdexcept>
#include <utility>
#include <tuple>

// template library headers
#include <Eigen/Eigen>

// independent headers
#include "xyz.h"
#include "io.h"

//!
//! moltool main function
//!
int main(int argc, char** argv)
{
try
{   
	io::arguments arg(argc,argv);
	io::switcharg<false> help(arg,"-h","--help", "print help message");
	io::switcharg<false> version(arg,"-v","--version", "print version info");
	io::vararg<std::string> selected_range (arg,"-s","--select","select atoms def.:'all'","all");
	io::vararg<std::string> iselected_range (arg,"-is","--inv-select","inverse select atoms def.:'none'","none");
	io::switcharg<false> sortz(arg,"-soz","--sort-z", "sorting in z direction");
	io::vararg<std::string> sortframesby(arg,"-sortby","--sortby","sort by variable <varname>","");
	io::switcharg<false> sortframesbynoa(arg,"-sortbynoa","--sortbynoa","sort frames by number of atoms");
	io::vararg<std::string> measure_dist(arg,"-d","--distance","distance of i,j atoms","");
	io::vararg<std::string> measure_zdiff(arg,"-dz","--diff-z","z coord difference of i,j atoms","");
	io::vararg<std::string> measure_angle(arg,"-angle","--angle","angle of i,j,k atoms","");
	io::vararg<std::string> measure_anglexy(arg,"-anglexy","--angle-xy-plane","measure angle of i,j atoms and the xy plane","");
	io::vararg<std::string> measure_dihedral(arg,"-dihedral","--dihedral","dihedral of i,j,k,l atoms","");
	io::vararg<std::string> measure_inversion(arg,"-inversion","--inversion","inversion of i,j,k,l atoms","");
	io::vararg<std::string> remove_var(arg,"-rv","--remove","remove <variable>","");
	io::vararg<std::string> match_value(arg,"-mv","--match-value","match value from a <list file>","");
	io::vararg<std::string> repeat(arg,"-r","--repeat","repeat unit cells in lattice vector direction <n1,n2,n3>","");
	io::vararg<std::string> nearest_atom(arg,"-na","--nearest-atom","n-th nearest selected atom to i, <n,i>","");
	io::switcharg<false> num_of_frames(arg,"-nof","--num-of-frames", "print number of frames"); 
	io::vararg<std::string> shake_atoms(arg,"-shake","--shake","shake selected atoms <n,dr>","");  //! ?
	io::vararg<std::string> stretch(arg,"-stretch","--stretch","stretch selected atoms with unitcell in z <minz,dz,maxz>","");
	io::multivararg<std::string> move(arg,"-move","--move","(mult) move selected atoms  <from,to,n> (from:x/y/z or i)");  
	io::multivararg<std::string> rotate(arg,"-rotate","--rotate","rotate selected atoms <o,a,angle,n> "); 
	io::vararg<std::string> alignz(arg,"-alignz","--alignz","align selected atoms into z direction by i,j atoms",""); 
	io::switcharg<false> split (arg,"-split","--split", "split frames"); 
	io::vararg<std::string> add_left(arg,"-al","--add-left","add <xyz> file to left side","");
	io::vararg<std::string> add_right(arg,"-ar","--add-right","add <xyz> file to right side","");
	io::vararg<int> snapshots(arg,"-snap","--snapshots","select about <n> frames",0);
	io::multivararg<std::string> write_variables(arg,"-wv","--write-var","(mult) write variables <varname>");
	io::switcharg<false> print_switch(arg,"-p","--print", "print frames"); 
        io::vararg<std::string> match_uc(arg,"-muc","--match-unitcell","match unitcell vector file","");

	io::posarg<std::string> input_fn(arg,"input file name", "");
	io::posarg<std::string> output_fn(arg,"output file name", "");
	
	if(*help)
	{
		arg.print_help(std::cout,"moltool");
		return EXIT_SUCCESS;
	}
	if(*version)
	{
		std::cout << "version 0.4" << std::endl;
		return EXIT_SUCCESS;
	}

	// input output streams
	io::input_stream input(*input_fn);
	io::output_stream output(*output_fn);
	
	std::istream &is = input.stream();
	std::ostream &os = output.stream();

	// load all xyz structure (100 atom with coord,veloc,force and 100000 frames ~ 1Gb in memory)
	std::vector<io::xyz<Eigen::Vector3d> > frames;
	io::xyz<Eigen::Vector3d> frame;
	while(is >> frame) frames.push_back(frame);

	// select range
	if((*selected_range == "all")&&(*iselected_range == "none"))
	{
		for(size_t i = 0; i < frames.size(); i++)
		for(int index = 1; index <= (int)frames[i].length(); index++)
		frames[i][index].mark = 1;
	}
	else
	{
		io::range_expression<int> re;
		if(*selected_range != "all")
		if(match("%1",*selected_range,&re))
		for(size_t i = 0; i < frames.size(); i++)
		frames[i].mark(*re);

		if(*iselected_range != "none")
		if(match("%1",*iselected_range,&re))
		{
			for(size_t i = 0; i < frames.size(); i++)
			for(int index = 1; index <= (int)frames[i].length(); index++)
			if(std::find((*re).begin(),(*re).end(),index) == (*re).end())
			frames[i][index].mark = 1; else frames[i][index].mark = 0;
		}
	}

	///// IN PLACE DATA ACTIONS ///////

	if(~sortframesby)  //! works for only double
	{
		for(size_t k = 0; k < frames.size(); k++)
		{
			size_t minl = k;
			for(size_t l = k+1; l < frames.size(); l++)
			if(frames[l].var_double(*sortframesby) < frames[minl].var_double(*sortframesby) )  minl = l;
			std::swap(frames[minl],frames[k]);
		}
	} else if(~sortframesbynoa) 
	{
		for(size_t k = 0; k < frames.size(); k++)
		{
			size_t minl = k;
			for(size_t l = k+1; l < frames.size(); l++)
			if(frames[l].length() < frames[minl].length() )  minl = l;
			std::swap(frames[minl],frames[k]);
		}
	} else if(*sortz)  
	{
		for(size_t k = 0; k < frames.size(); k++)
		{
			io::xyz<Eigen::Vector3d> & fr = frames[k];
			for(int i = 1; i <= fr.length(); i++)
			if(fr[i].mark == 1)
			{
				size_t minj = i;
				for(int j = i+1; j <= fr.length(); j++)
				if(fr[j].mark == 1)
				if(fr[j] < fr[minj]) minj = j;
				if(minj != i)
				fr.swap(minj,i);
			}
		}
	} else if(~measure_dist)
	{
		int i1,i2;
		if(io::match("%1,%2",',',*measure_dist,&i1,&i2))
		for(size_t i = 0; i < frames.size(); i++)
		{
			double d = (frames[i][i1].coord - frames[i][i2].coord).norm();
			std::ostringstream var_name;
			var_name << "d_" << i1 << "_" << i2;
			frames[i].var_double(var_name.str()) = d;
		}
	} else if(~measure_zdiff)
	{
		int i1,i2;
		if(io::match("%1,%2",',',*measure_zdiff,&i1,&i2))
		for(size_t i = 0; i < frames.size(); i++)
		{
			double d = (frames[i][i1].coord - frames[i][i2].coord).z();
			std::ostringstream var_name;
			var_name << "dz_" << i1 << "_" << i2;
			frames[i].var_double(var_name.str()) = d;
		}
	} else if(~measure_angle)
	{
		int i1,i2,i3;
		if(io::match("%1,%2,%3",',',*measure_angle,&i1,&i2,&i3))
		for(size_t i = 0; i < frames.size(); i++)
		{
			Eigen::Vector3d a  = frames[i][i1].coord-frames[i][i3].coord;
			Eigen::Vector3d b  = frames[i][i3].coord-frames[i][i2].coord;
			double d = acos(a.normalized().dot(b.normalized()));
			std::ostringstream var_name;
			var_name << "angle_" << i1 << "_" << i2 << "_" << i3;
			frames[i].var_double(var_name.str()) = d;
		}
	} else if(~measure_dihedral)
	{
		int i1,i2,i3,i4;
		if(io::match("%1,%2,%3,%4",',',*measure_dihedral,&i1,&i2,&i3,&i4))
		for(size_t i = 0; i < frames.size(); i++)
		{
			Eigen::Vector3d a = frames[i][i1].coord - frames[i][i2].coord;
			Eigen::Vector3d b = frames[i][i2].coord - frames[i][i3].coord;
			Eigen::Vector3d c = frames[i][i3].coord - frames[i][i4].coord;
			double cosphi = a.cross(b).normalized().dot(b.cross(c).normalized());
			double d =acos(cosphi); ; //!TODO
			std::ostringstream var_name;
			var_name << "dihedral_" << i1 << "_" << i2 << "_" << i3 << "_" << i4;
			frames[i].var_double(var_name.str()) = d;
		}
	} else if(~measure_inversion)
	{
		int i1,i2,i3,i4;
		if(io::match("%1,%2,%3,%4",',',*measure_inversion,&i1,&i2,&i3,&i4))
		for(size_t i = 0; i < frames.size(); i++)
		{
			Eigen::Vector3d a = frames[i][i1].coord - frames[i][i2].coord;
			Eigen::Vector3d b = frames[i][i1].coord - frames[i][i3].coord;
			Eigen::Vector3d c = frames[i][i1].coord - frames[i][i4].coord;

			double cosphi = a.cross(b).normalized().dot(c.normalized());

			double d = 0.5*2.1415 - acos(cosphi); //!TODO
			std::ostringstream var_name;
			var_name << "inversion_" << i1 << "_" << i2 << "_" << i3 << "_" << i4;
			frames[i].var_double(var_name.str()) = d;
		}
	} else if(~remove_var)
	{
		for(size_t i = 0; i < frames.size(); i++)
		frames[i].remove_var(*remove_var);
	} else if(~match_value)
	{
		std::vector<std::string> listlines;
		std::vector<double> vals;
		io::load(*match_value,"#",listlines);
		double p;
		for(std::vector<std::string>::iterator it = listlines.begin(); it != listlines.end(); it++)
		if(io::match("%1",*it,&p)) vals.push_back(p);
		if(frames.size() != vals.size())
			throw std::overflow_error("number of frames differs from value list size");
		for(size_t i = 0; i < frames.size(); i++)
			frames[i].var_double(*match_value)= vals[i];
	} else if(~match_uc)
	{
		std::vector<std::string> listlines;
		std::vector<Eigen::Vector3d> v1s;
                std::vector<Eigen::Vector3d> v2s;
                std::vector<Eigen::Vector3d> v3s;
		io::load(*match_uc,"#",listlines);
		double p1,p2,p3,p4,p5,p6,p7,p8,p9;
		for(std::vector<std::string>::iterator it = listlines.begin(); it != listlines.end(); it++)
		if(io::match("%1 %2 %3  %4 %5 %6  %7  %8  %9",*it,&p1,&p2,&p3,&p4,&p5,&p6,&p7,&p8,&p9))
                {
                    Eigen::Vector3d vtemp; 
                    
                    vtemp << p1,p2,p3;
                    v1s.push_back(vtemp);
                    
                    vtemp << p4,p5,p6;
                    v2s.push_back(vtemp);
                    
                    vtemp << p7,p8,p9;
                    v3s.push_back(vtemp);
                    
                }
		if(frames.size() != v1s.size())
			throw std::overflow_error("number of frames differs from uc list size");
		for(size_t i = 0; i < frames.size(); i++)
                {
			frames[i].a1() = v1s[i];
                        frames[i].a2() = v2s[i];
                        frames[i].a3() = v3s[i];
                }
	} else if(~nearest_atom)
	{
		int nth;
		int index;
		if(io::match("%1,%2",',',*nearest_atom,&nth,&index))
		for(size_t i = 0; i < frames.size(); i++)
		{
			std::set<int> nearest_set;
			int nearest;
			for(int n = 1;n <= nth; n++)
			{
				nearest = 0;
				for(int k = 1; k <= frames[i].length(); k++)
				if(frames[i][k].mark == 1)
				if(nearest_set.find(k) == nearest_set.end())
				{
					if(nearest == 0) nearest = k;
					if((frames[i][index].coord - frames[i][k].coord).norm() < 
						(frames[i][index].coord - frames[i][nearest].coord).norm()) nearest = k;
				}
				nearest_set.insert(nearest);
			}
			if(nearest == 0)
				throw std::logic_error("no n-th nearest out of k if k < n");
			
			std::ostringstream name;
			name << "%nearest_" << nth << "_" << index;
			frames[i].var_int(name.str()) = nearest; 
		}
	}

	//// OUTPUT ACTIONS ////

	// number of frames 
	if(*num_of_frames)
	{
		os << frames.size() << std::endl;
	} 
	else if(~shake_atoms)
	{
		int N;
		double dr;
		std::vector<io::xyz<Eigen::Vector3d> > outframes;
		if(io::match("%1,%2",',',*shake_atoms,&N,&dr))
		for(size_t i = 0; i < frames.size(); i++)
		{
			io::xyz<Eigen::Vector3d> fr = frames[i];
			outframes.push_back(fr);
			Eigen::Vector3d drv;
			for(int n = 0; n < N; n++)
			{
				for(int index = 1; index <= (int)fr.length(); index++)
				if(fr[index].mark == 1)
				{
					double num = static_cast<double>( rand() ) * (1.0-0.0) / static_cast<double>( RAND_MAX ) + 0.0; 
					drv.setRandom();
					drv.normalize();
					fr[index].coord = frames[i][index].coord + dr*num*drv; //!TODO
				}
				outframes.push_back(fr);
			}
		}
		
		for(size_t i = 0; i < outframes.size(); i++)
		os << outframes[i];
	} 
	else if(~stretch)
	{
		double dz;
		double mindz;
		double maxdz;
		std::vector<io::xyz<Eigen::Vector3d> > outframes;
		if(io::match("%1,%2,%3",',',*stretch,&mindz,&dz,&maxdz))
		for(size_t i = 0; i < frames.size(); i++)
		{
			io::xyz<Eigen::Vector3d> fr = frames[i];
			Eigen::Vector3d drv;
			double zz = mindz;
			while(zz <= maxdz)
			{
				drv << 0.0,0.0,zz;
				for(int index = 1; index <= (int)fr.length(); index++)
				if(fr[index].mark == 1)
					fr[index].coord = frames[i][index].coord + drv;
				fr.a3() = frames[i].a3() + drv;
				outframes.push_back(fr);
				zz+=dz;
			}
		}
		
		for(size_t i = 0; i < outframes.size(); i++)
		os << outframes[i];
	}
	else if(~repeat)
	{
		std::vector<io::xyz<Eigen::Vector3d> > outframes;

		int n1,n2,n3;
		if(io::match("%1,%2,%3",',',(*repeat),&n1,&n2,&n3)) ;
		else throw "syntax error in --repeat expression";

		for(size_t i = 0; i < frames.size(); i++)
		{
			const io::xyz<Eigen::Vector3d> & fr = frames[i];
			io::xyz<Eigen::Vector3d> fro = frames[i];
			
			Eigen::Vector3d a1 = fr.a1();
			Eigen::Vector3d a2 = fr.a2();
			Eigen::Vector3d a3 = fr.a3();


			for(int i1 = 0; i1<=n1;i1++) 
			for(int i2 = 0; i2<=n2;i2++)
			for(int i3 = 0; i3<=n3;i3++)
			if(i1+i2+i3 != 0)
			{
				for(int index = 1; index <= (int)fr.length(); index++)
				if(fr[index].mark == 1)
				{
					io::xyz<Eigen::Vector3d>::record_type atom = fr[index];
					atom.coord = (i1*a1+i2*a2+i3*a3)+atom.coord;
					fro += atom;
				}
			}

			fro.a1() = n1*a1;
			fro.a2() = n2*a2;
			fro.a3() = n3*a3;
			outframes.push_back(fro);

		}

		for(size_t i = 0; i < outframes.size(); i++)
		os << outframes[i];
		
	} 
	else if(~move)
	{
		std::vector<io::xyz<Eigen::Vector3d> > outframes;
		for(size_t n = 0; n < (*move).size(); n++)
		{
			std::string sb;
			std::string se;
			int N;
			if(io::match("%1,%2,%3",',',(*move)[n],&sb,&se,&N)) ;
			else if(io::match("%1,%2",',',(*move)[n],&sb,&se)) N = 0;
			else throw "syntax error in -move expression";

			Eigen::Vector3d vb,ve,u;
			double x,y,z;
			int index1 = 0;
			if(io::match("%1/%2/%3",'/',sb,&x,&y,&z)) vb << x,y,z;
			else if (io::match("%1",sb,&index1)) ;
			else throw "syntax error in -move expression";

			int index2 = 0;
			if(io::match("%1/%2/%3",'/',se,&x,&y,&z)) ve << x,y,z;
			else if (io::match("%1",se,&index2)) ;
			else throw "syntax error in -move expression";

			for(size_t i = 0; i < frames.size(); i++)
			{
				io::xyz<Eigen::Vector3d> fr = frames[i];
				
				if(index1 != 0) vb = fr[index1].coord;
				if(index2 != 0) ve = fr[index2].coord;

				u = ve-vb;
				u.normalize();
				double l = (ve-vb).norm();
				double d = l / (N > 1 ? N-1 : 0.5);

				Eigen::Vector3d drv;
				double ds = 0.0;
				if(N==0) ds = l;
				while(ds <= l)
				{
					drv = ds*u;
					for(int index = 1; index <= (int)fr.length(); index++)
					if(fr[index].mark == 1)
						fr[index].coord = frames[i][index].coord + drv;
					outframes.push_back(fr);
					ds += d;
				}
			}

		}
		
		for(size_t i = 0; i < outframes.size(); i++)
		os << outframes[i];
		
	} 
/*	else if(~alignz)
	{
		std::vector<io::xyz<Eigen::Vector3d> > outframes;
		for(size_t n = 0; n < (*move).size(); n++)
		{
			std::string sb;
			std::string se;
			if(io::match("%1,%2",',',(*move)[n],&sb,&se)) ;
			else throw "syntax error in -alignz expression";

			Eigen::Vector3d vb,ve,u;
			double x,y,z;
			int index1 = 0;
			if(io::match("%1/%2/%3",'/',sb,&x,&y,&z)) vb << x,y,z;
			else if (io::match("%1",sb,&index1)) ;
			else throw "syntax error in -move expression";

			int index2 = 0;
			if(io::match("%1/%2/%3",'/',se,&x,&y,&z)) ve << x,y,z;
			else if (io::match("%1",se,&index2)) ;
			else throw "syntax error in -move expression";

			for(size_t i = 0; i < frames.size(); i++)
			{
				io::xyz<Eigen::Vector3d> fr = frames[i];
				
				if(index1 != 0) vb = fr[index1].coord;
				if(index2 != 0) ve = fr[index2].coord;

				u = ve-vb;
				u.normalize();
				double l = (ve-vb).norm();
				double d = l / (N > 1 ? N-1 : 0.5);

				Eigen::Vector3d drv;
				double ds = 0.0;
				if(N==0) ds = l;
				while(ds <= l)
				{
					drv = ds*u;
					for(int index = 1; index <= (int)fr.length(); index++)
					if(fr[index].mark == 1)
						fr[index].coord = frames[i][index].coord + drv;
					outframes.push_back(fr);
					ds += d;
				}
			}

		}
		
		for(size_t i = 0; i < outframes.size(); i++)
		os << outframes[i];
		
	} */
	else if(~rotate)
	{
		std::vector<io::xyz<Eigen::Vector3d> > outframes;
		for(size_t n = 0; n < (*rotate).size(); n++)
		{
			std::string sax0;
			std::string sax1;
			int N;
			double a;

			if(io::match("%1,%2,%3,%4",',',(*move)[n],&sax0,&sax1,&a,&N)) ;
			else throw "syntax error in -rotate expression";

			Eigen::Vector3d vax0,vax1,ax;

			double x,y,z;
			int index1 = 0;
			if(io::match("%1/%2/%3",'/',sax0,&x,&y,&z)) vax0 << x,y,z;
			else if (io::match("%1",sax0,&index1)) ;
			else throw "syntax error in -rotate expression";

			int index2 = 0;
			if(io::match("%1/%2/%3",'/',sax1,&x,&y,&z))	vax1 << x,y,z;
			else if (io::match("%1",sax1,&index2)) ;
			else throw "syntax error in -rotate expression";

			for(size_t i = 0; i < frames.size(); i++)
			{
				io::xyz<Eigen::Vector3d> fr = frames[i];
				
				if(index1 != 0) vax0 = fr[index1].coord;
				if(index2 != 0) vax1 = fr[index2].coord;

				ax = vax1-vax0;
				ax.normalize();
				double da = a / (N > 1 ? N-1 : 0.5);

				Eigen::Transform3d rot;
				double an = 0.0;
				while(an <= a)
				{
					
					// TODO
					//rot = ;

					for(int index = 1; index <= (int)fr.length(); index++)
					if(fr[index].mark == 1)
						fr[index].coord = rot*(frames[i][index].coord-vax0)+vax0;
					outframes.push_back(fr);


					an += da;
				}
			}

		}
		
		for(size_t i = 0; i < outframes.size(); i++)
		os << outframes[i];
		throw "it is a beta version and this feature has not been implemented yet"; //!TODO
		
	}
	else if(~split)
	{
		for(size_t i = 0; i < frames.size(); i++)
		{
			std::ostringstream name;
			name << std::right << std::setfill('0') << std::setw(4) << i << ".xyz";
			std::ofstream ofst(name.str().c_str());
			ofst << frames[i];
			ofst.close();
		}
	} 
	else if(~add_left)
	{
		std::ifstream isleft((*add_left).c_str());

		io::xyz<Eigen::Vector3d> le;

		isleft >> le;
		Eigen::Vector3d lercorn = le[1].coord;
		for(int i = 1; i <= le.length(); i++)
		{
			if(le[i].coord.x() < lercorn.x() ) lercorn.x() = le[i].coord.x();
			if(le[i].coord.y() < lercorn.y() ) lercorn.y() = le[i].coord.y();
			if(le[i].coord.z() < lercorn.z() ) lercorn.z() = le[i].coord.z();
		}
		lercorn+=le.a3();

		isleft.close();

		std::vector<io::xyz<Eigen::Vector3d> > outframes;
		for(size_t k = 0; k < frames.size(); k++)
		{
			io::xyz<Eigen::Vector3d> fr = frames[k];
			
			Eigen::Vector3d lcorn = fr[1].coord;
			Eigen::Vector3d rcorn = fr[1].coord;
			for(int i = 1; i <= fr.length(); i++)
			{
				if(fr[i].coord.x() < lcorn.x() ) lcorn.x() = fr[i].coord.x();
				if(fr[i].coord.y() < lcorn.y() ) lcorn.y() = fr[i].coord.y();
				if(fr[i].coord.z() < lcorn.z() ) lcorn.z() = fr[i].coord.z();
			}
			//rcorn += fr.a3;
			
			io::xyz<Eigen::Vector3d> l;
			l = le;
			for(int i = 1; i <= l.length(); i++)
				l[i].coord += lcorn - lercorn;
			fr += l;
			fr.a3() = fr.a3() + le.a3();
			outframes.push_back(fr);
		}

		for(size_t i = 0; i < outframes.size(); i++)
		os << outframes[i];
	} 
	else if(~add_right)
	{
		std::ifstream isright(	(*add_right).c_str() );

		io::xyz<Eigen::Vector3d> ri;

		isright >> ri;
		Eigen::Vector3d rilcorn = ri[1].coord;
		for(int i = 1; i <= ri.length(); i++)
		{
			if(ri[i].coord.x() < rilcorn.x() ) rilcorn.x() = ri[i].coord.x();
			if(ri[i].coord.y() < rilcorn.y() ) rilcorn.y() = ri[i].coord.y();
			if(ri[i].coord.z() < rilcorn.z() ) rilcorn.z() = ri[i].coord.z();
		}
	
		isright.close();
		std::vector<io::xyz<Eigen::Vector3d> > outframes;
		for(size_t k = 0; k < frames.size(); k++)
		{
			io::xyz<Eigen::Vector3d> fr = frames[k];
			
			Eigen::Vector3d rcorn = fr[1].coord;
			for(int i = 1; i <= fr.length(); i++)
			{
				if(fr[i].coord.x() < rcorn.x() ) rcorn.x() = fr[i].coord.x();
				if(fr[i].coord.y() < rcorn.y() ) rcorn.y() = fr[i].coord.y();
				if(fr[i].coord.z() < rcorn.z() ) rcorn.z() = fr[i].coord.z(); 
			}
			rcorn += fr.a3();
			
			io::xyz<Eigen::Vector3d> r;
			r = ri;
			for(int i = 1; i <= r.length(); i++)
				r[i].coord += rcorn - rilcorn;
			fr += r;
			fr.a3() = fr.a3() + ri.a3();
			outframes.push_back(fr);
		}

		for(size_t i = 0; i < outframes.size(); i++)
		os << outframes[i];
	}
	else if(~snapshots)
	{
		int N = *snapshots;
		int every = ((N > 1) ? (frames.size() / (N-1)) : (frames.size()*2));
		every++;
		std::vector<io::xyz<Eigen::Vector3d> > outframes;
		for(size_t i = 0; i < frames.size(); i+=every)
			outframes.push_back(frames[i]);
		for(size_t i = 0; i < outframes.size(); i++)
		os << outframes[i];
	}
	else if(~write_variables)
	{
		os << '#';
		for(std::vector<std::string>::const_iterator it = (*write_variables).begin(); it != (*write_variables).end(); it++)
			os << std::left << std::setw(15) << std::setprecision(8) << std::fixed << *it << ' ';
		os << std::endl;
		for(size_t i = 0; i < frames.size(); i++)
		{
			for(std::vector<std::string>::const_iterator it = (*write_variables).begin(); it != (*write_variables).end(); it++)
			if((*it)[0]=='%')
				os << std::left << std::setw(15) << std::setprecision(8) << std::fixed << frames[i].var_int(*it) << ' ';
			else if ((*it)[0]=='$')
				os << std::left << std::setw(15) << std::setprecision(8) << std::fixed << frames[i].var_string(*it) << ' ';
			else
				os << std::left << std::setw(15) << std::setprecision(8) << std::fixed << frames[i].var_double(*it) << ' ';
			os << std::endl;
		}
	}
	else
	{
		if(*print_switch)
		{
			for(size_t i = 0; i < frames.size(); i++)
			print(std::cout,frames[i]);
			if( !output.is_std() )
			{
				for(size_t i = 0; i < frames.size(); i++)
				os << frames[i];
			}
		}
		else for(size_t i = 0; i < frames.size(); i++) os << frames[i];
	}
} catch(const char* message)
{
	std::cerr << message << std::endl;
	std::cerr << "(exit)" << std::endl;
	return EXIT_FAILURE;
} 
catch(const std::exception & e)
{
	std::cerr << e.what() << std::endl;
	std::cerr << "(exit)" << std::endl;
	return EXIT_FAILURE;
}
	return EXIT_SUCCESS;
}

